#include "initial_state.hpp"
#include "channels.hpp"
#include "hamilton.hpp"
#include <boost/numeric/ublas/io.hpp>

double f(int jj,int dx,int N,double Omega){
  double pi=3.141592653589793238462;
  return exp(-2*pi*double(jj)*dx/double(N))*sin(2*pi*double(jj)*Omega/double(N));
}

int main(){
  
  int N=40;
  int dx=1;
  double Omega1=9;
  double Omega2=7;
  double Omega3=2;
  
  vector<double> y1(N);
  vector<double> y2(N);
  vector<double> y3(N);
    
  for(int ii=0;ii<N;ii++){
    
    y1[ii]=f(ii,dx,N,Omega1);
    y2[ii]=f(ii,dx,N,Omega2);
    y3[ii]=f(ii,dx,N,Omega3);
  }
  rate r1(y1,double(dx));
  rate r2(y2,double(dx));
  rate r3(y3,double(dx));
  
  rate r4=r1;

  matrix<double,row_major,std::vector<double> > m1(2,2,{0,1,0,0});
  matrix<double,row_major,std::vector<double> > m2(2,2,{0,0,1,0});
  matrix<double,row_major,std::vector<double> > m3(2,2,{1,0,0,-1});
  
  channels chans;
  //Key 0 implicitly
  chans.set_op(m1);
  //Key 1 explicitlty
  chans.set_op(m2,1);
  //Key 2 implicitly
  chans.set_op(m3);
  //Key 0 implicitly
  chans.set_decay_rate(r1);
  //Key 1 explicitlty
  chans.set_decay_rate(r2,1);
  //Key 2 implicitly
  chans.set_decay_rate(r3);

  //Key 0 implicitly
  chans.set_lamb_rate(r1);
  //Key 1 explicitlty
  chans.set_lamb_rate(r2,1);
  //Key 2 implicitly
  chans.set_lamb_rate(r3);

  vector<std::complex<double> > a(4);
  initial_state is;
  initial_state is1(a);
  initial_state is2(is1);
  initial_state is3=is2;
  std::cout<<"Initial state  = "<<is()<<std::endl;
  std::cout<<"Initial state 1 = "<<is1()<<std::endl;
  std::cout<<"Initial state 2 = "<<is2()<<std::endl;
  std::cout<<"Initial state 3 = "<<is3()<<std::endl;
  
  //Hamiltonian
  hamilton H;
  hamilton H1(m3);
  hamilton H2=H1;
  hamilton H3(H1);

  std::cout<<"Number of channels: "<<chans.get_number_of_channels()<<"."<<std::endl;
  std::cout<<"Operator 0: "<<chans.get_op(0)<<"."<<std::endl;
  std::cout<<"Operator 0 times Rate 0 at time 4.12: "<<chans(0,4.12)<<"."<<std::endl;
  std::cout<<"Operator 1 "<<chans.get_op(1)<<"."<<std::endl;
  std::cout<<"Operator 1 times Rate 1 at time 4.12: "<<chans(1,4.12)<<"."<<std::endl;
  std::cout<<"Operator 2: "<<chans.get_op(0)<<"."<<std::endl;
  std::cout<<"Operator 2 times Rate 2 at time 4.12: "<<chans(2,4.12)<<"."<<std::endl;
  std::cout<<"\\sum_{k=1}^3 S_k(4.12) C_k^\\dag C_k: "<<chans(4.12)<<std::endl;
  std::cout<<"H = "<<H()<<std::endl;
  std::cout<<"H1 = "<<H1()<<std::endl;
  std::cout<<"H2 = "<<H2()<<std::endl;
  std::cout<<"H3 = "<<H3()<<std::endl;
  return 0;
}
